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Abstract. We believe the following three ingredients are enough to explain 
the mystery of the arrow of time: (1). equations of dynamics of gas molecules, 
(2). chaotic instabilities of the equations of dynamics, (3). unavoidable per- 
turbations to the gas molecules. The level of physical rigor or mathematical 
rigor that can be reached for such a theory is unclear. 



1. The Theory 

There have been a lot of writings on the arrow of time; nevertheless, the problem 
still remains a mystery. In the equilibrium thermodynamics context, the arrow of 
time refers to the second law of thermodynamics. Consider the example of gas in a 
box, each gas molecule moves according to the equation of dynamics (free flying and 
collision) which is reversible in the coordinate time; on the other hand, the whole 
body of gas can only move in an entropy increasing direction among thermodynamic 
equilibria, e.g. moving from half box to whole box. The entropy increasing direction 
can be viewed as an intrinsic time direction — the arrow of time. In fact, there are 
many similar irreversible features in nature; for instances, living things can only 
grow older, sugar and water do not spontaneously un-mix, broken glass does not 
spontaneously fix itself. Explaining all these arrows of time is a challenging task. 
According to L. Boltzmann [l]:"The second law (of thermodynamics) can never 
be proved mathematically by means of the equations of dynamics alone." Here we 
shall approach the problem with the following three ingredients; 

(1) equations of dynamics of gas molecules, 

(2) chaotic instabilities of the equations of dynamics, 

(3) unavoidable perturbations to the gas molecules. 

Trying to use chaos to explain irreversibility is not a new attempt, see for 
example I. Prigogine's writings [7] [9] |10) [8]. Of course, there are many other 
attempts from different perspectives. The most recent one is that of Maccone [4] 
from quantum perspective, see also its critiques [2] [3] [6] [5]. 

As Boltzmann said, the equations of dynamics alone will not be able to explain 
irreversibility. Unavoidable perturbations to the gas molecules have to be taken 
into account. These perturbations are inside the isolated system, and they can 
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be for examples physical noises and chemical fluctuations. If the dynamics of the 
gas molecules is non-chaotic, then these small perturbations will not generate any 
substantial effect; and the dynamics of the whole gas body will be reversible at 
least partially. The reality is that the system of equations of gas molecules is a 
non-integrable system and its dynamics is chaotic. Under such chaotic instabilities, 
the small unavoidable perturbations totally alter and liberate the dynamics of the 
whole gas body during the relaxation time from one thermodynamic equilibrium to 
another. In fact, the final thermodynamic equilibrium is independent of the initial 
conditions of the gas molecules, under such chaotic instabilities and unavoidable 
perturbations. For the example of gas in half box moving to full box, even the 
coordinate time and velocity of every gas molecule are reversed at any time during 
the relaxation time from half box to full box, the final thermodynamic equilibrium 
will still be the full box. Thus we arrive at the diagram for the arrow of time, 
Figure [T] 




Figure 1. The diagram of the arrow of time. The 'Past', 'Now' 
and 'Future' are coordinate time, and the 'Entropy' is the thermo- 
dynamic equilibrium entropy. 

Starting from the coordinate time 'Now', in both forward and backward co- 
ordinate time, the thermodynamic equilibrium entropy is increasing. The entropy 
time is folded into only positive direction — the arrow of time. In a more general- 
ized sense, for a living thing at 'Now', even its every molecule's coordinate time is 
reversed, the living thing will not get younger, rather continue to grow older! 

Next we will explain a little more on chaotic instabilities and the effect of un- 
avoidable perturbations. The main character of chaotic instabilities is 'sensitive 
dependence on initial data'. This is commonly characterized by amplifications of 
initial perturbations. For the traditional lower dimensional chaos, such amplifica- 
tions are exponential in time. For the dynamics of gas moleculesm such amplifica- 
tions are actually linear in time as shown later in the numerical simulation of one 
dimensional gas dynamics. The dynamics of the gas molecules is chaotic in both 
forward and backward coordinate time. Due to the unavoidable perturbations, the 
real physical orbits of the gas molecules are not the mathematical orbits prescribed 
by the equations of dynamics and the initial conditions. Only between two consec- 
utive perturbations among all gas molecules, the gas molecules follow a segment of 
the mathematical orbits. The duration between two consecutive perturbations is 
extremely small. So the real physical orbits quickly diverge from the mathematical 
ones. For an A'^ particles system on the interval [0, 1] for example, the phase space 
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has 2N dimensions, coordinated by the particles' locations and velocities. The re- 
gion where all particles are in the half interval [0, 1/2], has measure 2^^. The main 
theme of our study is to investigate the possibility of return to this small region 
when the time and velocities are reversed at some moment on the orbit originating 
from a point in this small region. 

2. A Popular Challenger to Irreversibility — The Poincare Recurrence 

Theorem 

The unwarranted argument here is that due to the Poincare recurrence theorem, 
full box gas will return near half box gas after enough time. Physicists making such 
an argument apply the Poincare recurrence theorem without its condition. The 
crucial condition is that one can find a bounded invariant region in the phase space 
of all the gas molecules coordinated by their locations and momenta. Then inside 
the bounded invariant region, Poincare recurrence theorem holds. The system of 
equations of gas molecules is a non-integrable system for which such a bounded 
invariant region does not exist. Even for near-integrable Hamiltonian systems, the 
well-known mechanism of Arnold diffusion can lead to unbounded drifting of the 
momenta. 

3. Ergodicity 

There are two types of chaotic dynamics of gas molecules when the whole 
system moves from one thermodynamic equilibrium to another. Type I chaotic 
dynamics is the 'transitional chaotic dynamics' during the relaxation time, while 
type II chaotic dynamics is the 'saturated chaotic dynamics' when the system is in 
a thermodynamic equilibrium. The transitional chaotic dynamics has a non-zero 
mean motion, and it is the non-zero mean motion that eliminates the possibility of 
ergodicity. The saturated chaotic dynamics has no mean motion and may have er- 
godicity. Since the system is in a thermodynamic equilibrium, its entropy maintains 
constant and is compatible with ergodicity of the gas molecules' dynamics. 

4. Numerical Simulations 

4.1. 1-Dimensional Dynamics. The main setup is as follows: Particles 
move in the interval [0,1], particles collide with each other and the walls (the 
two boundaries and 1 of the interval [0, 1]) elastically. In particular, after they 
collide, two particles exchange their speeds to maintain the conservation of momen- 
tum assuming all particles have the same mass; on the other hand, after it collides 
with a wall, a particle just reverses its velocity. 

4.1.1. Reversibility. In Figure [2j initially all the particles are randomly po- 
sitioned in the half interval [0,0.5], and the initial velocities of the particles are 
randomly chosen from [—0.5,0.5]. After each collision, a random perturbation of 
magnitude less than 10~^ is introduced to the velocities, while no perturbation is 
introduced to the positions of the particles. Special care needs to be taken when a 
speed is small after a collision. As time increases, the particles will evolve into the 
full interval [0, 1]. After 10,000 collisions, the particles velocities are reversed and 
the time is reversed too. Then after each new collision, check if the particles are all 
inside the half interval [0,0.5]. The return probability is given by the quotient 
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Figure 2. Return probability for particles initially randomly po- 
sitioned in the half interval [0, 0.5]. 

of the total number of collisions N,. after which the particles are all inside the half 
interval [0, 0.5] and 10, 000: 

Pr = A^,./10000 . 

For any fixed number of particles Np^ 10 such numerical experiments are conducted, 
and the average return probability (Pr) is ploted in Figure [2] 

In Figure [3j initially all the particles are randomly positioned in the full interval 
[0, 1], and the initial velocities of the particles are randomly chosen from [—0.5, 0.5]. 
No random perturbation is introduced. After each collision, check if the particles 
are all inside the half interval [0,0.5]. The numerical experiment ends after 10,000 
collisions. The return probability Pr is defined the same as above. The average 
return probability (Pr) of 10 such numerical experiments is ploted in Figure [Sj 

In Figure [4j initially all the particles are evenly (with equal spacing) positioned 
in the full interval [0,1], and the initial velocities of the particles are randomly 
chosen from [—0.5,0.5]. After each collision, a random perturbation of magnitude 
less than 10~^ is introduced to the velocities, while no perturbation is introduced 
to the positions of the particles. After each collision, check if the particles are 
all inside the half interval [0,0.5]. The numerical experiment ends after 50,000 
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Figure 3. Return probability for particles initially randomly po- 
sitioned in the full interval [0, 1]. 

collisions. The return probability Pr is defined the same as above (with 10, 000 
replaced by 50,000). The average return probability (Pr) of 10 such numerical 
experiments is ploted in Figure |4] 

From our preliminary numerical experiments, we conclude that 

• Without perturbation, no matter how many particles are initially po- 
sitioned in the half interval, and no matter how many total collisions, 
reversibility back to the half interval is always achieved as long as the 
computer accuracy is sufficient. This is due to the reversibility of the 
mathematical orbits. 

• With perturbation, reversibility back to the half interval is possible for 
small number of particles {Np < 16). When the number of particle Np > 
16, extremely tiny perturbations (of magnitude less than 10~^) will lead 
to irreversibility. This is due to the sensitive dependence on initial data 
of the orbits in the phase space. Take Np = 16 as the example, in this 
case, the phase space has 32 dimensions coordinated by the locations and 
velocities of the 16 particles. When the 16 particles are positioned in 
the half interval [0,0.5], and their initial velocities are chosen, a point 
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Figure 4. Return probability for particles initially evenly posi- 
tioned in the full interval [0, 1]. 



in the 32-diniensional phase space is chosen, from which a mathematical 
orbit exists. This rigorous mathematical orbit is reversible. But such 
a mathematical orbit has sensitive dependence on initial data. Even a 
tiny perturbation can lead to another diverging orbit. This causes the 
irreversibility of the numerical orbit when tiny random perturbations are 
introduced. We believe that the specifics of the perturbations and the 
specific ways in which they are introduced are important to irreversibility. 
What is crucial for irreversibility is the existence of perturbation. 

In Figure [5] initially all the particles are randomly positioned in the half in- 
terval [0,0.5], and the initial velocities of the particles are randomly chosen from 
[—0.5, 0.5]. After each collision, a random perturbation of magnitude less than 10~^ 
is introduced to the velocities, while no perturbation is introduced to the positions 
of the particles. Special care needs to be taken when a speed is small after a colli- 
sion. As time increases, the particles will evolve into the full interval [0, 1]. After 
different relaxation time (10000, 20000, 30000), the particles velocities are reversed 
and the time is reversed too. Then the particles are tracked the same amount of 
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Figure 5. Return (reverse) probability for different relaxation time. 



time in the reverse time, at the end check if the particles are all inside the half inter- 
val [0,0.5]. The new return probability Prv (or more precisely reverse probability) 
is defined by the quotient of the number of experiments with successful reverse N^y 
and the total number of experiments N: 

Prv = Nrv/N . 

Due to the size of the computation, we usually run the number of experiments N 
around 50 to a few hundreds. 

4.1.2. Chaotic Dynamics. First we look at the time series plots for two particles 
case (TVp = 2). In Figure [6j for four different initial conditions, the location time 
series of the right particle are plotted. These time series look chaotic. But such 
chaotic dynamics may not be in the chaos in the classical sense. Bear in mind that 
at a collision, the velocity has a discontinuity. Next we explore further the chaotic 
dynamics along the line of Liapunov exponents. 

If we introduce small perturbations only on the initial locations of the particles, 
this will cause the perturbed collisions and the unperturbed collisions happen at 
a little different time, then the velocity deviation can be order one during a small 
time interval inbetween the perturbed and the unperturbed collisions. On the other 
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Figure 6. The location time series of the right particle in the two 
particles system. 

hand, such a velocity deviation does not contribute substantially to the location 
deviation. So the location deviation maintains the same order as the initial loca- 
tion perturbation. Our numerical experiments verified this. Therefore, the main 
contribution for irreversibility comes from velocity perturbations. 

If we introduce small perturbations only on the initial velocities (or both on 
the initial velocities and on the initial locations) of the particles, then the initial 
velocity deviation will generates a location deviation which is linear in time, with 
speed given by the initial velocity deviation. Near collisions, the location deviation 
will diverge from this main line, resulting in a rough oscillating continous non- 
differentiable location deviation function in time. The velocity deviation mainly 
follows the horizontal line in time given by the initial velocity deviation. Near 
collisions, it will just to other horizontal lines of order one away from the horizontal 
axis. Overall, the velocity deviation is a piecewise continuous function in time, 
with discontinuities at the the collisions, see Figures [7|8j For Np = 2, there are 
5 horizontal lines for velocity deviation (Figure l7|. The main line is the one very 



ON THE ARROW OF TIME 



9 



2.Q 



40 



60 



ao 



100 



Figure 7. Velocity deviation in Euclidean norm. Np = 2. The 
larger dots mark the velocity deviation, and the smaller dots mark 
the exact horizontal lines. Initial perturbations are on both veloc- 
ity and position. Each perturbation is bounded by 10^^. 



close to the horizontal line, given by the initial velocity deviation of order 10 ^: 

where (ui, V2) and {v[,v'2) are the unperturbed and perturbed velocities of the two 
particles. The other four lines are near: 

2vi, 2V2, \/2|vi+U2|, V2\vi-V2\. 

For Np = 3, there are 10 horizontal lines for velocity deviation. For Np — 4, there 
are 17 horizontal lines for velocity deviation. The location deviations are shown 
in Figures 9]|T0 The slope of the main line in the location deviation is given by 



the main line of the velocity deviation, which is very close to the horizontal line. 
The scattered points below the the main line in the location deviation come from 
the other lines of the velocity deviation, in lieu of collisions. Notice that, all the 
scattered points are below the main line in the location deviation. This can be 
explained from a closer look at the collision. Let vi and V2 be the velocities of 
the unperturbed particles, vi + -Oi and V2 + V2 be the velocities of the perturbed 
particles. Without loss of generality, we assume that the unperturbed particles 
collide earlier. Then there are essentially two possibilities as shown in Figure [TT] 
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Figure 8. Velocity deviation in Euclidean norm. iVp = 16. Initial 
perturbations are on both velocity and position. Each perturbation 
is bounded by 10^^. 

In the case of Figure [Tl|a), right after the collision, both particle 1 and both particle 
2's location deviations will decrease. In the case of Figure [Tl|b), right after the 
collision, particle 1 (the circle) 's location deviations will increase, while particle 2 
(the square) 's location deviations will decrease. More specifically, right after the 
collision, particle I's velocity deviation is 

V2 - {v\ + i)i), 
while particle 2's velocity deviation is 

V\ - (W2 + ^2)- 

Therefore, particle I's velocity has an extra deviation (due to the collision) of 
while particle 2's velocity has an extra deviation of 
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Figure 9. Location deviation in Euclidean norm. Np = 2. Initial 
perturbations are on both velocity and position. Each perturbation 
is bounded by 10^^. 



Thus the extra gain (due to the collision) of particle I's location deviation is the 
same with the extra loss of particle 2's location deviation. In term of ^i-norm 
(sum of every particle's location deviation), in the case of Figure [njb), the total 
location deviation of particles 1 and 2 is unchanged in lieu of the collision. On the 
other hand, in term of £2-norm (Euclidean norm), the total location deviation may 
increase or decrease. So far, we never observed the 'increase' case. If a particle 
collides with a boundary, it is easy to see that the particle's location deviation will 
decrease right after the earlier collision. 

Scattered points below the main line in the location deviation always appear 



when the experiment time is long enough, see Figure 12 for the two particles case. 
If we increase the size of the initial perturbation, scattered points appear earlier. In 
Figures 13 and 14 we use the initial perturbation size of 10~^ and the experiment 
time of 5000 and 10000, we can see a global picture in time of the location deviation 
for the two particles case. Bear in mind that the blue region represents a continuous 
non-differentiable function in time. 

Ill conclusion, what we have observed is that the particle dynamics has sensitive 
dependence upon initial data, but instead of "exponential" sensitive dependence 
like for traditional chaos, we observed here "linear" sensitive dependence; that is, 
the particle location deviation induced by an initial perturbation grows linearly 
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Figure 10. Location deviation in Euclidean norm. Np — 16. 
Initial perturbations are on both velocity and position. Each per- 
turbation is bounded by 10~^. 
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Figure 11. Essential types of collisions. 



in time, instead of exponentially in time. It is such a linear growth that leads to 
irreversibility. 
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Figure 12. Location deviation in Euclidean norm. Np = 2. Ini- 
tial perturbations are on both velocity and position. Each pertur- 
bation is bounded by 10^^. 



4.2. 2-Dimensional Dynamics. The main setup is now as follows: Circular 
disks of radius r (instead of points in ID case) move in the rectangle (of size 
8x4), the elastic collision principle is applied to the disks colliding with each other 
and with the walls of the rectangle. The locations and velocities of the disks are 
identified by those of their centers. No rotational effect is considered. 

4.2.1. Reversibility. Limited numerical experiments indicates that the return 
probability Pr (defined above) from the half to the whole rectangle follows the law 
of Pr ~ 2~^ where N is the number of disks; while in the ID case, the return 
probability does not exactly follow this law due to particles blocking each other. 

4.2.2. Chaotic Dynamics. In the 2D case, the dynamics of the disks is com- 
pletely chaotic. 

First (Figure 15) we choose the number of disks = 6, and the radius of each 
disk r = 0.05. The initial velocities of the disks are randomly chosen with speed 
bounded by 10. The initial locations of the disks are also randomly chosen within 
the rectangle. By adding a perturbation of size 10~^^ to the initial locations, we 
observe the location deviations of the disks in time. Figure 15 shows the time series 
of the location deviations in Euclidean norm. 

Next we compute the Liapunov exponent. For a perturbation of size 10~^^ 
to the initial locations, we stop the computation when the Euclidean norm of the 
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Figure 13. Location deviation in Euclidean norm. Np = 2. Ini- 
tial perturbations are on both velocity and position. Each pertur- 
bation is bounded by 10^"*. 



location deviations reaches 10 ^ . The time step of the computation is 0.002. For 
each of the last 100 time steps, we compute the quantity 

At Ap{t) ' 

where Ap represents the Euclidean norm of the location deviations, and At rep- 
resents the time step. Then we average these 100 quantities to get an estimated 
Liapunov exponent for this numerical experiment. Finally, we average these esti- 
mated Liapunov exponents for 200 such numerical experiments to obtain the final 
estimated Liapunov exponent. The dependence of the final estimated Liapunov 



exponent upon the size and the number of the disks is shown in Figure 16 



Finally the slide shows on the dynamics of six disks are shown in Figures 17 
andfTl 



5. Epilogue 

We are convinced that the three ingredients mentioned in the section on Theory 
are enough to explain the mystery of the arrow of time. The level of physical rigor 
or mathematical rigor that can be reached for this theory is unclear. The numerical 
simulations above (especially in the 2D case) offer positive support to our Theory. 
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Figure 14. Location deviation in Euclidean norm. Np = 2. Ini- 
tial perturbations are on both velocity and position. Each pertur- 
bation is bounded by 10^"*. 



A huge amount of writings on the arrow of time exists in the literature, and we 
made no attempt to scan through them. 
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Figure 15. Time scries of the location deviations in Euclidean 
norm in the 2D case. 
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Figure 16. The final estimated Liapunov exponent as a function 
of the size and the number of the disks. 
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Figure 17. The evolution of the six disks and the evolution of the 
perturbed six disks. Since the perturbation size is 10~^, initially 
the unperturbed and the perturbed disks coincide almost com- 
pletely. The radius of the disk is 0.25, and the rectangle domain is 
8x4. 
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(d) t = 5.248 



Figure 18. The same setup as in Figure [T7| except that the per- 
turbation size is 10^^^ 



